% function to detrend, de-bias, calibrate, fft, and cross & auto power calculation
% References Bendat & Piersol, 1971, Metronix 1997
% included sub functions parzen.m and doparzen.m in this file
% to save computing time 30.4.3
% checked, checked & checked the parzening pl do not do it again
% tried to bypass parzen averaging with simple median - very bad results 22.5.3
% Robust band average good results 5.6.3
% latest date 16.6.3

function[SPMatrix,ProcDef] = GetSpectra(TS,ProcDef),

sensor = ProcDef.sensor;
Notchi = TS.filter.NotchF(2:5,:);		% NOTE ! there may be problems 
IsNotch = TS.filter.Filt(2:5,:);		%        if notch do not have right dimension..
oi = size(Notchi);
if oi(2) == 1,				% Nasty way. But it just works
	notch = Notchi(find(IsNotch)); 
elseif oi(2) == 2,
 	[x,y]=find(IsNotch);
   notch = Notchi(1:oi(2),1:oi(2));
end;
hpb1 = TS.filter.Filt(1);			% NOTE ! high pass filter is always disabled
prob = TS.probe.dist;	% 	 in calib05.m
gain = TS.probe.gain;
nband = TS.head.band;
nchannel = TS.head.nchannel;
block = ProcDef.block;
TLine1 = ProcDef.TLine1;
TLRad1 = ProcDef.TLRad1;
TLFreq1 = ProcDef.TLFreq1;
nfrq=ProcDef.nfrq;
notchdef = ProcDef.notchdef;
Array = ProcDef.Array;

[TIndex,len] = GetSelArray(ProcDef);
ProcDef.len=len;
ProcDef.SelStacks = ones([1,len]);
ProcDef.Weights = ones([nfrq,len]);

if len < 1,
   fprintf('At least one stack required - error in GetSpectra\n');
   return;
end;


disp('GetSpectra -> Calibration...');

for i = 1:nchannel,     %|
   if length(notch) < 1,%|
      notch = [0]; 		%| inserted 25.3.3
   end;						%| reason - in some rts fls notch = [0 0 0 0]
   
   nt=find(notch(:,notchdef(i)));
   if length(nt)>0,
      nt=notch(nt);
   end;
   [cplx,frq1] = calib05(sensor(i),nband,block,prob(i),nt,hpb1,gain(i));
     
CalData(i,:) = cplx;
end;

%------------native declarations---------------
disp('GetSpectra -> Now making matrices of cross and auto spectra....');

WI = hanning(block);


SPMatrix = zeros([nfrq len nchannel nchannel]);

fff = len/20;
aaa = 1;
fprintf('0                 100percent\n');
for i = 1:len,
   
   if i/fff > aaa,
      fprintf('#');
      aaa = aaa+1;
   end;
   
   
% FFT, calibration etc.

for u = 1:nchannel,
   K = CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),u),0);
	Channel1(:,u) = fft(K.*WI)/block;
   Channel(:,Array(u)) = Channel1(2:(block/2)+1,u).*conj(CalData(u,1:block/2)'); %
end; 	
% OOOPS the ' transpose made all the probelm. Solved it with another conj -30,10,02 (b'day eve)
% auto and cross spectra---------------------------
for ii = 1:nfrq,
   par = parzen1(TLRad1(ii)); % Bit inefficient (as it is called many times to do the same job)
   SPMatrix(ii,i,:,:) = CompuSPM1(Channel,par,TLine1(ii),TLRad1(ii));
end;
%------------------------------------------------
%clear Channel;

end;
fprintf('\n Over..Switching to main\n');

%helper functions------------------------------

function[SPM] = CompuSPM1(Channel, par,TL,TR),

a = size(Channel);
nchannel = a(2);

for i = 1:nchannel,
   for j = 1:i,
%talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i))); % just summing 
talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i)).*par); % summing w.r.t parzen window%
%talkal = rbpz(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i))); % Robust Band averaging
      if i ~= j,
       SPM(i,j) = real(talkal);
       SPM(j,i) = (-1)*imag(talkal); 
      elseif i == j,
       SPM(i,j) = talkal;
     end;
  end;
end;

function[PF] = parzen1(radius),
    if radius > 1,

%---- compute norm-factor for parzening 

      nFactor = 1.0; 
		PF(1) = 1.0;
      for line = 2:radius,

        u = (line/(radius+1))*pi;
        v = (sin(u)/u)^2.0;
        PF(line) = v;
        nFactor = nFactor+(2*v*v);
      end;
      nFactor =sqrt(nFactor);

%      -- parzening line 1 

      PF(1) = 1.0/nFactor;

%     -- parzening line 2 to radius      
	for line = 2:radius,
        PF(line) = PF(line)/nFactor;
     end;
     ParInt = sum(PF(2:radius))*2 + PF(1);
PF = PF/ParInt;
else,
   PF = 1;
   end;
   PF= [PF(radius:-1:2) PF]';
   
   
   %----trying to make parzen average a robust est 25.5.3''
   
function[C] = rbpz(S),

sig_m = 1.483*median(abs(S-median(S))); % median absolute deviation 
c_m = 1.5*sig_m;								% Scale
k = find(abs(S) <= c_m);					% assigning Huber weights 
w_m(k) = 1;
q = find(abs(S) > c_m);
w_m(q) = c_m./abs(S(q));					
C = sum(S.*w_m')/sum(w_m);  			% Weighted X spectra
